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1.  Introduction 


The  mechanism  and  rate  of  vibrational  energy  transfer  in  gas-surface  collisions  has 
been  of  great  interest  both  experimentally  and  theoretically.  Experimentally,12  it 
has  been  observed  that  accommodation  coefficients  for  high-frequency  molecular 
vibrations  tend  to  be  small,  whereas  corresponding  coefficients  for  molecules  with 
low-frequency  vibrations  close  to  the  impact  surface  are  comparatively  much 
higher  (although  some  exceptions  have  been  reported3).  Theoretically,  gas-surface 
interactions  have  been  studied  using  a  variety  of  simulation  techniques.  Alexander 
et  al.4  employed  a  scattering  model  to  study  energy  transfer  for  atoms  impacting 
self-assembled  monolayer  surfaces.  Lucchese  and  Tully5  employed  a  stochastic 
trajectory  approach  to  study  the  impact  of  nitric  oxide  molecules  on  the  001  surface 
of  the  lithium  fluoride  crystal  as  a  function  of  vibrational  frequency,  surface 
temperature,  and  translational  energy.  In  that  work,  it  was  found  that  the  dominant 
mechanism  of  energy  transfer  between  the  surface  and  impacting  gas  was  coupling 
of  vibrations  to  surface  phonons  and  that  a  relatively  small  amount  of  energy 
resulted  from  repartitioning  of  the  vibrational  energy  in  rotational  and  translation 
modes. 

In  this  work,  using  atomistic  molecular  dynamics  (MD)  simulation,  the  mechanism 
of  energy  deposition  by  a  shocked  diatomic  gas  (nitrogen  [N2])  into  a  stationary 
target  is  studied  as  a  function  of  multiple  variables  including  gas  density,  impact 
velocity,  and  target  rigidity.  We  focus  here  primarily  on  the  resulting  gas  dynamics 
and  detail  the  partitioning  of  the  energy  among  the  available  rotational  and 
vibrational  channels  as  a  function  of  impact  condition.  The  results  suggest  that 
rotational  excitations  are  important  at  all  impact  velocities  and  that  vibrational 
excitation  in  the  gas  is  only  important  for  high-velocity  impacts.  The  results  also 
suggest  that  the  rate  of  energy  deposition  into  the  vibrational  channels  of  the  gas  is 
a  function  of  the  density. 

2.  Computational  Methods 


2.1  Simulation  Cell 

A  general  schematic  of  the  simulation  cells  used  in  this  work  is  shown  in  Fig.  1. 
Each  cell  consists  of  a  target  placed  on  the  left  side  of  the  cell  and  a  diatomic  gas 
(N2)  occupying  the  right  side  of  the  cell  where  the  Z  Cartesian  axis  is  normal  to  the 
target  surface  as  shown.  The  cell  dimensions  for  all  simulations  were  100.8  x  100.8 
x  820  A. 
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Fig.  1  Schematic  of  simulation  cell  (image  not  drawn  to  scale) 

For  each  simulation,  the  target  consisted  of  point  masses  (14.0  g/mol)  with  an  inter- 

o 

particle  separation  of  1.2  A.  The  target  layers  were  regularly  spaced  from  Z  =  0  to 

o 

Z  =  14.4  A  yielding  91,728  particles  in  the  target  region  of  the  simulation  cell.  Each 
pair  of  atoms  in  the  target  was  connected  by  springs  to  all  its  nearest  neighbors  and 
interacted  via  a  harmonic  potential 

V(r)taraet  =\k(r-r0)2,  (1) 

where  r  is  the  inter-particle  separation  (bonded  interactions  only)  and  r<>  is  the 

o 

equilibrium  bond  length.  In  this  work,  the  parameter  r0  equaled  1.2  A,  which  is 
equivalent  to  the  value  of  the  lattice  spacing  used  to  construct  the  target  layers. 
Therefore,  the  particles  in  the  target  already  occupied  their  equilibrium  lattice 
positions,  by  construction,  and  as  a  result,  there  was  minimal  reorientation  of  the 
particles  in  the  target  layer  during  cell  equilibration.  The  uniformity  of  the  target 
layer  geometry,  irrespective  of  target  rigidity  (to  be  discussed),  is  important  since 
this  removes  a  potential  source  of  variation  in  the  observed  energy  deposition  and 
gas  dynamics  as  the  impact  conditions  are  varied. 

The  force  constant  k  in  Eq.  1  controls  the  strength  of  the  interaction  between  bonded 
pairs  and  consequently  determines  the  target  stiffness.  Six  targets  were  prepared 
with  force  constant  values  of  100,  500,  1500,  3000,  6000,  and  10,000 
(kcal  mol-1  A-2)  and  as  a  rough  estimate  of  the  surface  stiffness  for  each  k,  one  can 
compute  the  vibrational  frequency  between  a  pair  of  target  particles  acting  as 
harmonic  oscillators  that  are  decoupled  from  their  other  neighbors.  In  this 
approximation,  the  frequency,  v,  between  each  target  atom  pair  is  given  by 


where  p  is  the  reduced  mass  between  the  particles  and  k  is  the  force  constant  value 
for  the  target.  The  resulting  frequencies  for  the  6  force  constants  used  in  this  work 
are  given  in  Table  1.  The  frequencies  range  from  a  low  energy  mode  of  410  cm-1 
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(characteristic  of  a  relatively  soft  target)  up  to  a  maximum  of  4100  cm  a  high- 
energy  vibration  characteristic  of  a  rigid  target. 

Table  1  Vibrational  frequencies,  v,  for  each  force  constant  k.  Frequencies  in  cm-1,  force 
constants  in  kcal  moF1  A-2. 


k 

V 

(cm-1) 

too 

410 

500 

917 

1500 

1588 

3000 

2246 

6000 

3176 

10000 

4100 

Following  placement  of  the  target  in  the  simulation  cell,  the  impacting  gas  was 
generated  by  population  of  the  empty  region  of  the  cell  above  the  target  with 
randomly  oriented  N2  molecules  at  multiple  experimental  densities,  p.  The 
experimental  density  of  nitrogen  gas  is  0.001251  g/cm3.6  This  is  a  factor  of  1000 
times  lower  than  the  densities  of  materials  (such  as  energetics  with  p  ~  1.8  g/cm3) 
that  are  usually  studied  in  atomistic  shock  simulations.  At  the  experimental  density 
and  cell  dimensions  used  in  this  work  (100.8  x  100.8  x  820  A),  only  215  N2 
molecules  are  present  in  the  gas  region  of  the  cell  (Fig.  2)  and  it  was  not  expected 
that  any  meaningful  physics  regarding  gas  dynamics  would  be  observed  with  this 
paucity  of  atoms. 


P  gas 


Fig.  2  Simulation  cells  with  various  N2  densities.  N2  atoms  are  blue  and  the  target  region  is 
red. 
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Due  to  the  low  particle  count  of  N2  at  the  experimental  gas  density,  simulation  cells 
with  N2  molecules  at  the  experimental  liquid  density  (puq  =  0.8084  g/cm3)  were 
prepared.  In  addition,  simulation  cells  using  values  of  0.1puq  and  0.5puq  were 
constructed  (Fig.  2)  to  study  the  effect  of  N2  density  on  energy  deposition  and 
redistribution.  Table  2  lists  the  total  simulation  size  (including  the  91,728  target 
atoms)  as  a  function  of  N2  density.  The  simulation  size  ranges  from  92,158  atoms 
for  the  simulation  cell  at  the  gas  density  (consisting  almost  entirely  of  target  atoms) 
up  to  369,776  atoms  for  the  simulation  cell  at  the  full  liquid  density.  Combining  the 
6  targets  with  the  4  N2  densities  yielded  24  different  simulation  cells  for  analysis. 

Table  2  Number  of  atoms  in  simulation  cell  (including  target)  as  function  of  N2  density 


p 

Number  of  atoms 

Pair 

92,158 

O.lpLiq 

119,380 

0.5pLiq 

230,752 

PLiq 

369,776 

2.2  Interaction  Potentials 

The  target-target  particle  interactions  were  governed  by  a  harmonic  potential 
(Eq.  1)  as  previously  described.  The  N2  molecules  were  treated  using  the 
Condensed-phase  Optimized  Molecular  Potentials  for  Atomistic  Simulation 
Studies  (COMPASS)7  force  field  where  the  intramolecular  interactions  of  N2  were 
given  by 

V(jyntra  =  c^r  _  boy  +  c^r  _  by3  +  cyr  _  (3) 

and  the  intermolecular  interactions  between  N2  molecules  were  of  the  Lennard- 
Jones  form 


V(r) 'X‘‘r  =  £  [2  (ff  -  3  (f)‘] .  (4) 

The  parameters  appearing  in  Eqs.  3  and  4  (c,  b0,  c,  o)  were  taken  from  the  literature.7 
The  functional  form  of  the  nitrogen-target  particle  interactions  can  be  chosen  freely, 
and  for  simplicity,  these  interactions  were  also  represented  by  Eq.  4  with  the  same 
parameter  values. 

2.3  Simulation  Protocol 

All  simulations  were  run  using  the  Large-scale  Atomic/Molecular  Massively 
Parallel  Simulator  (LAMMPS)8  software  package.  Each  simulation  cell  was 
equilibrated  in  the  isothermal-isochoric  (NVT)  ensemble  for  a  period  of  500  ps 
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using  a  1-fs  time  step  at  a  temperature  (T)  of  100  K.  The  simulation  employed 
periodic  boundary  conditions  in  the  X  and  Y  directions  (Fig.  1)  and  was 
nonperiodic  along  Z.  To  keep  the  target  region  from  translating  (a  necessity  for  the 
impact  studies),  the  bottom  layer  of  target  atoms  at  Z  =  0  was  kept  fixed  during  the 
simulations.  Diffusion  of  the  N2  molecules  from  the  “open”  end  of  the  simulation 
cell  (upper  edge  opposite  the  target)  was  prevented  by  using  the  “wall/reflect” 
option  in  LAMMPS.  With  this  technique,  if  a  nitrogen  atom  crossed  the  upper  cell 
boundary  by  a  distance  A,  it  was  placed  back  inside  the  cell  by  the  same  A  and  the 
sign  of  the  Z-component  of  the  velocity  for  that  particle  was  flipped,  thereby 
confining  all  atoms  to  the  interior  of  the  cell. 

After  equilibration,  N2  was  driven  into  the  target  at  multiple  impact  velocities,  Vimp. 
Momentum  of  N2  toward  the  target  was  introduced  by  addition  of  Vimp  to  the 
Z-component  of  the  velocity  vector,  vz,  for  each  nitrogen  atom.  The  values  of  vx 
and  vy  were  unchanged.  For  each  of  the  24  simulation  cells,  4  values  of  Vimp  (0.5, 
1.0,  1.5,  and  2.0  km/s)  were  used  yielding  96  simulations  total.  For  each  cell,  the 
impact  trajectory  was  integrated  in  the  microcanonical  (NVE)  ensemble  for  a  period 
of  100  ps.  However,  only  the  data  in  the  first  40  ps  were  used  for  determination  of 
total  energy  deposition  and  rovibrational  analysis  because  of  the  wall/reflect  option 
in  LAMMPS  as  previously  described.  This  option,  although  useful,  represents  a 
perturbation  on  the  translational  and  rovibrational  energies  since  it  alters  the 
molecular  geometry  and  velocities.  Observation  of  the  simulation  trajectories 
showed  that  after  driving  N2  into  the  target,  the  molecules  diffused  back  toward  the 
free  end  of  the  cell  and  reached  the  upper  edge  after  approximately  40  ps.  Therefore, 
analyses  of  the  gas  dynamics  were  terminated  at  that  point  so  that  the  results  would 
not  be  impacted  by  the  artificial  reflections  occurring  at  the  far  end  of  the  cell. 

2.4  Energy  Deposition  into  Target 

The  energy  of  the  target  was  computed  at  every  simulation  step  and  included  kinetic 
and  potential  energy  contributions  where  the  potential  energy  of  the  target  was 
given  by  summation  of  Eq.  1  over  all  bonded  pairs  in  the  target  region.  The  energy 
deposition,  Edep,  is  given  by 

Edep  —  ^Target, i  ^Target, o  >  (4) 

where  Erarget.iis  the  energy  of  the  compressed  target  at  the  /th  time  step  and  ETarget.o 
is  the  equilibrium  energy  of  the  target  before  impact.  The  target  energies,  ETarget,i 
were  computed  by  averaging  over  500  energy  values  centered  about  the  current 
time  step. 
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2.5  Energy  Distribution  in  N2 


The  partitioning  of  energy  in  N2  was  monitored  by  computing  the  kinetic  energy  of 
vibration  and  the  kinetic  energy  of  rotation  of  all  N2  molecules.  The  molecular 
rovibrational  energies  were  obtained  using  a  postprocessing  program,  developed  by 
the  author,  that  executes  the  following  steps  for  each  N2  molecule: 

1)  Compute  the  center  of  mass  and  translate  to  the  origin. 

2)  Subtract  velocity  contributions  corresponding  to  center  of  mass  translation. 

3)  Compute  the  inertia  tensor  and  diagonalize  yielding  the  principal  axes. 

4)  Rotate  the  local  molecular  coordinates  (and  atomic  velocities)  of  the  current 
molecule  to  the  principal  axis  orientation  using  the  convention  that  the 
eigenvector  with  a  zero  moment  of  inertia  (bond  axis  of  a  diatomic)  is 
oriented  along  the  X-axis. 

5)  Compute  the  kinetic  energy  of  vibration,  Kv  =  mv* . 

1  7  1  o 

6)  Compute  the  rotational  kinetic  energy,  Krot  =  -ly0Jy  +  -/zoj|,  where  Iy  and 
Iz  are  the  inertial  moments  (eigenvalues  obtained  in  step  3)  and  the  angular 

V 

velocities,  co,  about  each  axis  are  given  by  -  where  v  is  the  velocity  about 

the  axis  and  r  is  the  distance  of  the  atom  from  the  origin  in  the  local 
coordinate  system. 

3.  Results 


3.1  Energy  Deposition 

Representative  plots  of  the  change  in  target  energy  and  temperature,  relative  to  the 
equilibrium  state,  as  a  function  of  impact  velocity  for  the  target  with  k  =  1500  are 
presented  in  Figs.  3-6.  As  shown,  at  a  fixed  density,  the  energy  transferred  to  the 
target  increases  as  a  function  of  the  impact  velocity,  as  expected.  The  energy 
deposition  also  increases  as  a  function  of  the  density  of  N2  and  after  40  ps,  the 
amount  of  energy  transferred  to  the  target  for  the  0.5 puq  and  full  liquid  density  N2 
simulations  is  essentially  equivalent  (=20,000  kcal/mol).  This  suggests  saturation 
with  respect  to  density  for  these  small-length  scales.  The  temperature  rise  of  the 
target  ranges  from  approximately  4  K  for  impacts  of  N2  at  the  gas  density  to 
approximately  50  K  for  simulations  at  the  full  liquid  density.  However,  these 
energy  and  temperature  values  correspond  to  quantities  computed  40  ps  after 
impact  (the  point  at  which  the  simulations  were  terminated  due  to  the  issues 
discussed  previously)  and  the  final  energy  and  temperature  values  could  change  if 
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a  longer  simulation  with  a  larger  cell  was  used.  However,  based  on  the  slopes  of 
the  curves  after  40  ps  of  simulation  time,  it  is  not  expected  that  the  final  energy  and 
temperature  would  be  “significantly”  different  since  the  slopes  of  the  curves  are 
approaching  zero  in  most  cases. 
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Fig.  3  Change  in  target  energy  and  temperature,  relative  to  the  equilibrium  state,  as  a 
function  of  impact  velocity  (km/s)  for  k  =  1500  target  and  N2  at  the  experimental  gas  density 
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Delta  E  (kcal/mol) 


Fig.  4  Change  in  target  energy  and  temperature,  relative  to  the  equilibrium  state,  as  a 
function  of  impact  velocity  (km/s)  for  k  =  1500  target  and  O.lpuq  density  for  N2 
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Delta  E  (kcal/mol) 


Fig.  5  Change  in  target  energy  and  temperature,  relative  to  the  equilibrium  state,  as  a 
function  of  impact  velocity  (km/s)  for  k  =  1500  target  and  0.5puq  density  for  N2 
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Fig.  6  Change  in  target  energy  and  temperature,  relative  to  the  equilibrium  state,  as  a 
function  of  impact  velocity  (km/s)  for  k  =  1500  target  and  pLiq  density  for  N2 


Tables  3-8  report  the  energy  deposition  rate  (in  units  of  gigawatts  [GW])  per  mol) 
and  the  associated  temperature  increase  for  all  96  combinations  of  target,  N2 
density,  and  impact  velocity  studied  in  this  work.  Examination  of  the  data  shows 
that  the  energy  deposition  rate  or  “power”  (energy  per  unit  time)  follows  anticipated 
trends.  At  constant  density  and  constant  target  stiffness,  the  power  increases  as  a 
function  of  impact  velocity.  As  an  example,  for  the  k  =  1500  target  and  0.5puq  N2 
density,  the  power  increases  from  a  value  of  7.4  GW  for  Vi,np  =  0.5  km/s  to  a  value 
of  133.9  GW  for  Vimp  =  2.0  km/s.  Similarly,  at  constant  impact  velocity  and  constant 
density,  the  power  decreases  as  a  function  of  target  rigidity.  As  an  example,  at 
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Vimr  =  2.0  km/s  and  0.5puq  N2  density,  the  energy  deposition  rate  decreases  from 
384  GW  for  the  very  soft  target  (k  =  100)  to  only  61  GW  for  the  very  stiff  target 
(. k  =  10,000)  at  the  same  density.  Similar  trends  follow  from  other  target,  density, 
and  Vimp  combinations. 

Table  3  Energy  deposition  rate  into  target  with  k  =  100,  VimP(km/s),  AE  (kcal/mol),  AT  (K), 
and  rate  (GW/mol) 


k  =  100 


pair 

O.lpLiq 

Vimp 

AE 

AT 

Rate 

Vimp 

AE 

AT 

Rate 

0.5 

2340.4 

4.6 

14.0 

0.5 

2677.8 

5.3 

16.0 

1.0 

2742.9 

5.4 

16.4 

1.0 

8201.6 

16.3 

49.0 

1.5 

3306.8 

6.5 

19.8 

1.5 

19087.0 

38.0 

114.0 

2.0 

3809.6 

7.5 

22.8 

2.0 

36011.0 

71.6 

215.2 

0.5pLiq 

piiq 

0.5 

3318.1 

6.7 

19.8 

0.5 

3387.7 

6.8 

20.2 

1.0 

12876.6 

25.9 

76.9 

1.0 

13059.7 

26.3 

78.0 

1.5 

33322.2 

66.8 

199.1 

1.5 

32382.3 

64.9 

193.5 

2.0 

64264.0 

128.2 

384.0 

2.0 

62777.8  125.1 

375.1 

Table  4  Energy  deposition  rate  into  target  with  k  = 
and  rate  (GW/mol) 

k  =  500 

500,  Vimp  (km/s),  AE  (kcal/mol),  AT  (K), 

Pair 

O.lpLiq 

Vimp 

AE 

AT 

Rate 

Vimp 

AE 

AT 

Rate 

0.5 

1638.1 

3.2 

9.8 

0.5 

1541.6 

3.0 

9.2 

1.0 

2132.1 

4.2 

12.7 

1.0 

4617.1 

9.1 

27.6 

1.5 

2386.6 

4.7 

14.3 

1.5 

10543.1 

20.9 

63.0 

2.0 

3311.5 

6.6 

19.8 

2.0 

19887.7 

39.4 

118.8 

0.5pLiq 

piiq 

0.5 

1860.1 

3.6 

11.1 

0.5 

1772.2 

3.6 

10.6 

1.0 

6777.0 

13.4 

40.5 

1.0 

6632.9 

13.3 

39.6 

1.5 

18326.1 

36.4 

109.5 

1.5 

17308.6 

34.5 

103.4 

2.0 

37057.2 

73.5 

221.4 

2.0 

34683.3 

68.9 

207.2 
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Table  5  Energy  deposition  rate  into  targets  with  k  =  1500,  Vimp  (km/s),  AE  (keal/mol),  AT  (K). 
and  rate  (GW/mol) 


k  =  1500 


Pair 

O.lpLiq 

Vimp 

AE 

AT 

Rate 

Vimp 

AE 

AT 

Rate 

0.5 

1399.6 

2.8 

8.4 

0.5 

969.8 

1.9 

5.8 

1.0 

1668.6 

3.3 

10.0 

1.0 

2960.3 

5.9 

17.7 

1.5 

1899.8 

3.8 

11.4 

1.5 

6463.6 

12.8 

38.6 

2.0 

2404.5 

4.8 

14.4 

2.0 

12512.8 

24.8 

74.8 

0.5pLiq 

pliq 

0.5 

1235.0 

2.4 

7.4 

0.5 

1210.1 

2.4 

7.2 

1.0 

4311.4 

8.5 

25.8 

1.0 

4260.2 

8.5 

25.5 

1.5 

11233.0 

22.3 

67.1 

1.5 

10232.8 

20.4 

61.1 

2.0 

22404.4 

44.4 

133.9 

2.0 

20931.9 

41.6 

125.1 

Table  6  Energy  deposition  rate  into  targets  with  k  = 
and  rate  (GW/mol) 

k  =  3000 

3000,  Vimp  (km/s),  AE  (kcal/mol),  AT  (K), 

Pair 

O.lpLiq 

Vimp 

AE 

AT 

Rate 

Vimp 

AE 

AT 

Rate 

0.5 

877.4 

1.7 

5.2 

0.5 

861.8 

1.4 

5.1 

1.0 

922.0 

1.8 

5.5 

1.0 

2255.5 

4.2 

13.5 

1.5 

1338.5 

2.6 

8.0 

1.5 

4904.8 

9.4 

29.3 

2.0 

1760.6 

3.4 

10.5 

2.0 

9720.6 

18.9 

58.1 

0.5pLiq 

Pliq 

0.5 

777.6 

1.7 

4.6 

0.5 

786.2 

1.3 

4.7 

1.0 

3410.3 

6.9 

20.4 

1.0 

3407.9 

6.6 

20.4 

1.5 

8661.7 

17.3 

51.8 

1.5 

8141.2 

16.0 

48.6 

2.0 

17421.2 

34.6 

104.1 

2.0 

16278.1 

32.2 

97.3 

Table  7  Energy  deposition  rate  into  targets  with  k  = 
and  rate  (GW/mol) 

k  =  6000 

6000,  Vimp 

(km/s),  AE  (kcal/mol),  AT  (K), 

pair 

O.lpLiq 

Vimp 

AE 

AT 

Rate 

Vimp 

AE 

AT 

Rate 

0.5 

820.7 

1.6 

4.9 

0.5 

368.2 

0.7 

2.2 

1.0 

843.8 

1.7 

5.0 

1.0 

1208.3 

2.4 

7.2 

1.5 

1067.1 

2.1 

6.4 

1.5 

3191.1 

6.3 

19.1 

2.0 

1218.9 

2.4 

7.3 

2.0 

6633.8 

13.1 

39.6 

0.5pLiq 

Pliq 

0.5 

317.8 

0.7 

1.9 

0.5 

213.3 

0.4 

1.3 

1.0 

1936.5 

3.9 

11.6 

1.0 

1762.5 

3.5 

10.5 

1.5 

6364.2 

12.7 

38.0 

1.5 

6295.8 

12.5 

37.6 

2.0 

13559.3 

26.9 

81.0 

2.0 

12824.4 

25.4 

76.6 
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Table  8  Energy  deposition  rate  into  targets  with  k  =  10,000,  VimP  (km/s),  AE  (kcal/mol), 
AT  (K),  and  rate  (GW/mol) 


k  =  10,000 

pair 

O.lpLiq 

Vimp 

AE 

AT 

Rate 

Vimp 

AE 

AT 

Rate 

0.5 

457.6 

0.9 

2.7 

0.5 

150.3 

0.3 

0.9 

1.0 

627.7 

1.2 

3.8 

1.0 

553.7 

1.1 

3.3 

1.5 

669.8 

1.3 

4.0 

1.5 

1842.9 

3.6 

11.0 

2.0 

917.9 

1.8 

5.5 

2.0 

4222.5 

8.3 

25.2 

0.5pLiq 

P'i<i 

0.5 

92.9 

0.2 

0.6 

0.5 

39.5 

0.1 

0.2 

1.0 

828.4 

1.6 

4.9 

1.0 

690.0 

1.4 

4.1 

1.5 

3887.8 

7.7 

23.2 

1.5 

3697.6 

7.4 

22.1 

2.0 

10235.4 

20.3 

61.2 

2.0 

10328.3 

20.5 

61.7 

3.2  Energy  Redistribution  in  Gas 

Plots  of  the  molecular  vibrational  and  rotational  kinetic  energies,  summed  over  all 
No  molecules  in  the  cell,  for  the  0. 1  puq,  0.5puq,  and  puq  densities  ( k  =  1500  target) 
are  presented  in  Figs.  7-9.  These  plots  are  representative  of  the  results  obtained  for 
all  other  target,  density,  and  Vimp  combinations  except  those  for  the  experimental 
gas  density  that  were  indeterminate  due  to  the  low  number  of  N2  molecule.  As  such, 
this  experimental  density  is  not  discussed  further.  As  shown,  for  each  density,  the 
vibrational  kinetic  energy  remains  essentially  constant  (thermal  oscillations  aside) 
for  the  3  lowest  Vimp  values  suggesting  that  for  these  weaker  impacts,  there  is  no 
appreciable  change  in  vibrational  energy  within  the  40-ps  time  frame  of  these 
simulations.  However,  for  Vimp  =  2  km/s,  there  is  a  marked  increase  in  the  kinetic 
energy  of  vibration  suggesting  that  there  is  an  impact  threshold  upon  which 
vibrational  excitation  (either  via  direct  energy  deposition  or  through  rovibrational 
coupling)  becomes  important. 
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Fig.  7  Total  rotational  and  vibrational  kinetic  energies  for  k  =  1500  target  and  0.1  puq  N2 
density  as  a  function  of  impact  velocity  (km/s) 
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Rotational  Kinetic  Energy  (J/mol)  Vibrational  Kinetic  Energy  (J/mol) 
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Fig.  8  Total  rotational  and  vibrational  kinetic  energies  for  k  =  1500  target  and  0.5pLk|  N2 
density  as  a  function  of  impact  velocity  (km/s) 


Approved  for  public  release;  distribution  is  unlimited. 


15 


P|_iq 


o 


Time  (fs) 


Fig.  9  Total  rotational  and  vibrational  kinetic  energies  for  k  =  1500  target  and  puq  N2 
density  as  a  function  of  impact  velocity  (km/s) 

The  data  represented  in  Figs.  7-9  also  suggest  that  the  rate  of  energy  transfer  into 
the  vibrational  modes  is  a  function  of  the  density.  In  each  figure,  the  impact  portion 
of  the  simulation  begins  at  10  ps  with  the  steps  prior  to  that  representing  the  final 
steps  of  the  equilibration  phase.  For  the  0.1  pyq  density  simulation  (Fig.  7),  the 
vibrational  kinetic  energy  for  Vimp  =  2  km/s  begins  to  increase  at  approximately 
20  ps  (i.e.,  10  ps  after  the  initial  impact).  The  increase  occurs  at  approximately 
12  ps  for  the  0.5puq  simulation  (Fig.  8)  (i.e.,  2  ps  after  impact),  and  at  full  density 
(Fig.  9)  the  kinetic  energy  increase  begins  almost  instantaneously.  On  the  contrary, 
the  rotational  kinetic  energy  for  all  impact  and  density  values  shows  an 
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instantaneous  increase  immediately  after  target  impact  at  t  =  10  ps,  as  one  would 
expect,  since  rotation  is  a  much  lower  energy  process  than  vibration.  The  change  in 
rotational  kinetic  energy  for  Vimp  =  0.5  km/s  is  much  less  than  that  observed  for  the 
higher  impact  velocities  (masked  in  part  by  the  scale  of  the  plot).  This  is  due  to  the 
reduced  number  of  N2  collisions  with  the  target  resulting  from  this  comparatively 
weak  value  of  Vimp. 

Figures  10-12  are  renderings  of  the  simulation  cell  containing  the  k  =  1500  target 
(vimp  =1.5  km/s)  where  the  local  regions  of  the  cell  are  colored  by  the  average 
rovibrational  temperature  per  N2  molecule.  In  the  figures,  the  target  is  represented 
by  the  plane  of  atoms  at  the  far  left  edge  of  the  cell  and  each  colored  bin  is  of 

o 

dimension  5  x  5  x  15  A.  The  bins  in  the  target  region  remain  “cold”  in  the 
renderings  since  the  rovibrational  analysis  was  limited  to  the  N2  region  of  the  cell. 
The  times  used  to  label  the  snapshots  in  each  figure  are  relative  to  the  start  of  the 
impact.  At  t  =  0  (no  impact),  the  rovibrational  energy  is  equally  distributed 
throughout  N2  with  an  average  rovibrational  temperature  of  approximately  100  K 
per  local  bin.  However,  the  snapshots  taken  2.5  ps  after  impact  show  an  increase  in 
the  rovibrational  temperature  at  the  interface,  and  the  radial  extent  of  this  “hot  spot” 
appears  to  increase  as  a  function  of  the  N2  density.  At  t  =  25  ps,  the  gas  reaches  its 
maximum  compression  and  begins  to  flow  back  down  the  tube,  remaining 
rovibrationally  hot  through  t  =  40  ps  at  which  point  the  simulation  trajectory  was 
terminated.  Figures  10-12  also  suggest  that  the  internal  temperature  of  N2  after 
impact  is  a  function  of  density.  For  the  0.1  puq  density  simulation  (Fig.  10), 
molecular  internal  temperatures  exceed  1000  K  (red  bins)  after  impact.  However, 
when  the  density  is  increased  to  0.5puq  and  puq,  the  molecular  temperatures  only 
rise  to  approximately  500  K  and  300  K,  respectively.  This  is  likely  attributable  to 
an  increased  level  of  energy  dissipation  occurring  after  impact  between  the  N2 
molecules  that  results  from  the  increased  number  of  intermolecular  interactions  at 
the  higher  mass  densities.  This  may  also  be  why  the  size  of  the  hot  spot  increases 
with  density. 
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Fig.  10  Local  rovibrational  temperatures  for  O.lpuq  N2  density.  Times  measured  relative  to 
start  of  impact.  Target  is  at  far  left  of  cell. 
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Fig.  11  Local  rovibrational  temperatures  for  0.5puq  N2  density.  Times  measured  relative  to 
start  of  impact.  Target  is  at  far  left  of  cell. 


Approved  for  public  release;  distribution  is  unlimited. 

19 


0  ps 


2.5  ps 


I  1000K 
500K 


40  ps 


Fig.  12  Local  rovibrational  temperatures  for  puq  N2  density.  Times  measured  relative  to 
start  of  impact.  Target  is  at  far  left  of  cell. 

4.  Conclusion 


In  this  work,  the  physics  of  energy  deposition  and  rovibrational  excitation  was 
studied  using  MD  simulation.  The  results  suggest  that  rotational  excitations  are 
important  at  all  impact  velocities.  However,  vibrational  excitation  is  only  important 
for  high-velocity  impacts  (2  km/s  for  the  systems  in  this  work).  The  results  also 
suggest  that  the  rate  of  energy  deposition  into  the  vibrational  channels  is  a  function 
of  the  density.  The  simulations  conducted  in  this  work  used  nonreactive  potentials. 
However,  for  the  lower  density  simulations,  internal  particle  temperatures  exceeded 
1000  K,  which  in  principle  could  drive  molecular  dissociation  that  may  influence 
the  observed  mechanisms  of  energy  transfer.  Although  the  simulation  sizes  used  in 
this  work  are  beyond  the  application  of  quantum  mechanical  approaches,  the 
available  empirical  models  that  include  chemical  reactivity  (at  increased 
computational  expense,  however)  may  provide  insight  regarding  the  importance  of 
bond  dissociation  and  ionization  on  energy  deposition  and  partitioning. 
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